
clear
%%%%%%%%%%%%%%%%%%%%% Define parameters and starting values
cd 'C:\Users\at449\Desktop\replication package\matlab code'   %set working directory

vIncome=cell2mat(struct2cell(load('matrixIncomeG.mat'))); %load income


%default parameters
beta_I= -1.2 *10^(-7);
beta_r=5.8*10^(-3);

%loan demand parameters
phi=1.376e-05;
alphatilde=-1.679e-02;

vbetatilde= [1.486*10^(-03);-4.365e-06;1.500e-03]; %LTV, fees, teaser period

%loan choice parameters
sI=size(vIncome);
alpha_i=-4.5*10^(-5)*vIncome; % number of borrowers in row, rate in column

LTV_coef=2.3*10^(-1)+ normrnd(0,1*10^(-2),[sI(1),1]);
fees_coef= repmat(-9*10^(-4),sI(1),1);
teaser_coef=normrnd(0,2,[sI(1),1]);

%Fixed cost parameter
fixedcosts=1000;




beta_i=[LTV_coef,fees_coef,teaser_coef] ;  %  randn(10, 3); number of borrowers in row, nb of product characteristics in column

%marginal costs
mcLTV=-2.330e-01;
mcfees=-2.373e-04;
mcteaser=7.044e-02;
vmc=[mcLTV;mcfees;mcteaser];   % randn(3, 1);

mcLTV2=1.559e-03;  
mcfees2=-4.591e-08;  
mcteaser2=3.977e-03;
vmc2=[mcLTV2;mcfees2;mcteaser2];   %randn(3, 1);

bankmc=1.047e+01;   %mc

%%%%%%%%%%%%%%%%%%%%% start maximization problem

%initial products offered
mProd=cell2mat(struct2cell(load('matrixproducts.mat'))); %load product

mX_c=mProd(:,1:3)'; % one product per column
vr_c=mProd(:,4); % one product per line



%%%% check for other starting point
WTP=-LTV_coef./ alpha_i ; %just to check WTP for LTV

mProdbis=mProd; %check another initial point
mProdbis(mProdbis(:,1)<=90,1)=90; %update LTV
mProdbis(:,4)=mProdbis(:,4)+((mProdbis(:,1)-mProd(:,1))'*min(WTP))';% update rate check if this is smaller than 90 LTV rate for similar product

%initialize the bar u variable (i.e. market equilibrium)
Mutil=util(beta_i,alpha_i,mX_c,vr_c);
baru = PU(Mutil); % expected utility that each borrower get in the current data


%%tranform the initial product as a vector
x_0=mProd(:);
x_0bis=mProdbis(:);
smP=size(mProd);

A = [];
b = [];
Aeq = [];
beq = [];
lb = []; % may need to put a lower bound to rates
ub = []; %may need to put an upper bound to fees and LTV

%Solve max problem 1:define the variable that need to be optimized over: all varaible
%continuous
fun = @(x) objfn(x, vIncome, beta_I, beta_r, phi, alphatilde, vbetatilde,vmc,vmc2, bankmc,beta_i,alpha_i,smP);
nonlcon = @(x) constr(x,baru, beta_i,alpha_i,smP);
nonlcon2= @(x) constr2(x,baru, beta_i,alpha_i,smP);

mProd2=mProd;
mProd2(:,4)=mProd(:,4)-0.9;
x_0b=mProd2(:);
RSsolwc(:,4)=0;
x_0bb=RSsolwc(:);

[x3,fval3] = fmincon(fun,x_0,A,b,Aeq,beq,lb,ub,nonlcon2); 
[x3b,fval3b] = fmincon(fun,x_0bis,A,b,Aeq,beq,lb,ub,nonlcon2); 


%%%%%%Solve max problem 2:define the variable that need to be optimized over: fix the product then
%maximize rate

%create a list of potential product and solve for the optimal rate for each
%menus

ListProd={[90,0,5; 95,0,5 ;90,1000,5;95,1000,5;90,0,2; 95,0,2 ;90,1000,2;95,1000,2],[95,1000,5;90,0,2; 95,0,2 ;90,1000,2;95,1000,2]};
%prealocate list
ListRate={0,0};
ListProfits=[0,0];

for n = 1:2
%initialize the product (Can run the problem with different product and select the max)
mProdF= ListProd{n};
s=size(mProdF);
%initialize the price
x_0FC=zeros(s(1),1);

funFP = @(x) objfnFP(x, vIncome, beta_I, beta_r, phi, alphatilde, vbetatilde,vmc,vmc2, bankmc,beta_i,alpha_i,mProdF);
nonlconFP = @(x) constrFP(x,baru, beta_i,alpha_i,mProdF);
nonlcon2FP= @(x) constr2FP(x,baru, beta_i,alpha_i,mProdF);

%maximizaton problem
[xFP,fvalFP] = fmincon(funFP,x_0FC,A,b,Aeq,beq,lb,ub,nonlcon2FP); 
Profitsnet=fvalFP-fixedcosts*s(1);
ListRate{n}=xFP
ListProfits(n)=Profitsnet
end

%find the best menu
[maxValue, maxIndex] = max(ListProfits);
productmax=ListProd{maxIndex};
ratemax=ListRate{maxIndex};


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Save outcome into
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% dataset (model 2)
mX_cFC=productmax'; % one product per column
vr_cFC=ratemax;% prices

%%calculate the main variable as function of X, r

Mtheta=theta(vIncome, vr_cFC, beta_I, beta_r); % row contain borrower i column product c
MLoan=LoanFP(vIncome,phi, alphatilde, vbetatilde, mX_cFC,vr_cFC);% row contain borrower i column product c
vmc_c = MC(vmc,vmc2, bankmc, mX_cFC); % row contain  product c
si=size(beta_i); %to get the nb of borrowers si(1)
Mr=repmat(vr_cFC',si(1),1);                % interest rate matrix (for each i c )
Mmc=repmat(vmc_c,1,si(1));      % marginal costs matrix (for each i c )
MutilFC=util(beta_i,alpha_i,mX_cFC,vr_cFC) ; % utility if i chooses c
MSFC=marketshare(MutilFC); %generate the probability that i chooses c

nbpFC=sum(MSFC,1)';%number of people that chose a given product c
dataFC= [mX_cFC',vr_cFC,nbpFC];

% Specify the filename
filename = 'dataFC.csv';
% Save the matrix to CSV
writematrix(dataFC, filename);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Save outcome into
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% dataset (model 1)
%check solution
RSsol=reshape(x3,[smP(1),smP(2)]);
RSsol1=reshape(x3b,[smP(1),smP(2)]);
RSsol0_1=reshape(x_0bis,[smP(1),smP(2)]);


mean(RSsol0_1(:,1))
mean(RSsol1(:,1))
mean(RSsol0_1(:,4))
mean(RSsol1(:,4))


L=Loan(vIncome,phi, alphatilde, vbetatilde, mX_c,vr_c)

mProd2=reshape(x3b,[smP(1),smP(2)]);
mX_c=mProd2(:,1:3)'; % one product per column
vr_c=mProd2(:,4);

Mtheta=theta(vIncome, vr_c, beta_I, beta_r); % row contain borrower i column product c
MLoan=Loan(vIncome,phi, alphatilde, vbetatilde, mX_c,vr_c);% row contain borrower i column product c
vmc_c = MC(vmc,vmc2, bankmc, mX_c); % row contain  product c
si=size(beta_i); %to get the nb of borrowers si(1)
Mr=repmat(vr_c,1,si(1));                % interest rate matrix (for each i c )
Mmc=repmat(vmc_c,1,si(1));      % marginal costs matrix (for each i c )
Mutil=util(beta_i,alpha_i,mX_c,vr_c) ; % utility if i chooses c
MS=marketshare(Mutil); %generate the probability that i chooses c


nbp=sum(MS,1)';%number of people that chose a given product c
data= [mProd2,nbp];
data2=mProd;

% Specify the filename
filename = 'data.csv';
% Save the matrix to CSV
writematrix(data, filename);

% Specify the filename
filename = 'dataI.csv';
% Save the matrix to CSV
writematrix(data2, filename);



%%calculate the firm's profit

pi_ic=(Mtheta.* Mr- Mmc); % profits per loan units (i row c columns)
Pi=MS.*MLoan.*pi_ic ;% profits per loan units (i row c columns) * loan amount* pobability choose product
fval = sum( Pi , 'all' );% total profits

%Screening externality

mProd3=reshape(x3,[smP(1),smP(2)]);
mX_c3=mProd2(:,1:3)'; % one product per column
vr_c3=mProd2(:,4);


-((fval3-fval3b)/sI(1))./mean(L);
-((fval3-fval3b)/sum(mProd(:,1)<90))./mean(L);% 15 basis points
u1=util(beta_i,alpha_i,mX_c3,vr_c3);
u2=util(beta_i,alpha_i,mX_c,vr_c);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%% Define functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%generate firm's profits
function fval = objfn(x, vIncome, beta_I, beta_r, phi, alphatilde, vbetatilde,vmc,vmc2, bankmc,beta_i,alpha_i,smP)

%reshape x into product and rates
mProd2=reshape(x,[smP(1),smP(2)]);
mX_c=mProd2(:,1:3)'; % one product per column
vr_c=mProd2(:,4);% prices

%%calculate the main variable as function of X, r

Mtheta=theta(vIncome, vr_c, beta_I, beta_r); % row contain borrower i column product c
MLoan=Loan(vIncome,phi, alphatilde, vbetatilde, mX_c,vr_c);% row contain borrower i column product c
vmc_c = MC(vmc,vmc2, bankmc, mX_c); % row contain  product c
si=size(beta_i); %to get the nb of borrowers si(1)
Mr=repmat(vr_c,1,si(1));                % interest rate matrix (for each i c )
Mmc=repmat(vmc_c,1,si(1));      % marginal costs matrix (for each i c )
Mutil=util(beta_i,alpha_i,mX_c,vr_c) ; % utility if i chooses c
MS=marketshare(Mutil); %generate the probability that i chooses c


%%calculate the firm's profit

pi_ic=(Mtheta.* Mr- Mmc); % profits per loan units (i row c columns)
Pi=MS.*MLoan.*pi_ic ;% profits per loan units (i row c columns) * loan amount* pobability choose product
fval = -sum( Pi , 'all' );% total profits

end

% firms profits for fixed products
function fval = objfnFP(x, vIncome, beta_I, beta_r, phi, alphatilde, vbetatilde,vmc,vmc2, bankmc,beta_i,alpha_i,mProdF)

%reshape product
mX_c=mProdF'; % one product per column
vr_c=x;% prices

%%calculate the main variable as function of X, r

Mtheta=theta(vIncome, vr_c, beta_I, beta_r); % row contain borrower i column product c
MLoan=LoanFP(vIncome,phi, alphatilde, vbetatilde, mX_c,vr_c);% row contain borrower i column product c
vmc_c = MC(vmc,vmc2, bankmc, mX_c); % row contain  product c
si=size(beta_i); %to get the nb of borrowers si(1)
Mr=repmat(vr_c',si(1),1);                % interest rate matrix (for each i c )
Mmc=repmat(vmc_c,1,si(1));      % marginal costs matrix (for each i c )
Mutil=util(beta_i,alpha_i,mX_c,vr_c) ; % utility if i chooses c
MS=marketshare(Mutil); %generate the probability that i chooses c


%%calculate the firm's profit

pi_ic=(Mtheta.* Mr- Mmc'); % profits per loan units (i row c columns)
Pi=MS.*MLoan.*pi_ic ;% profits per loan units (i row c columns) * loan amount* pobability choose product
fval = -sum( Pi , 'all' );% total profits

end

%utility constraints
function [cons, conseq] = constr(x,baru, beta_i,alpha_i,smP)

%reshape x into product and rate
mProd2=reshape(x,[smP(1),smP(2)]);
mX_c=mProd2(:,1:3)'; % one charac per column
vr_c=mProd2(:,4);

Mutil=util(beta_i,alpha_i,mX_c,vr_c) ;% get utilities
pu = PU(Mutil);% get expected utility

hold_pc= pu-baru; %need to be a row vector

conseq = [hold_pc];

cons = []; % need to specify empty inequality constraints

end

%utility constraints fixed product
function [cons, conseq] = constrFP(x,baru, beta_i,alpha_i,mProdF)

%reshape product
mX_c=mProdF'; % one charac per column
vr_c=x;

Mutil=util(beta_i,alpha_i,mX_c,vr_c) ;% get utilities
pu = PU(Mutil);% get expected utility

hold_pc= pu-baru; %need to be a row vector

conseq = [hold_pc];

cons = []; % need to specify empty inequality constraints

end

%utility constraints
function [cons, conseq] = constr2(x,baru, beta_i,alpha_i,smP)

%reshape x into product and rate
mProd2=reshape(x,[smP(1),smP(2)]);
mX_c=mProd2(:,1:3)'; % one charac per column
vr_c=mProd2(:,4);

Mutil=util(beta_i,alpha_i,mX_c,vr_c) ;% get utilities
pu = PU(Mutil);% get expected utility

hold_pc= pu-baru; %need to be a row vector

conseq = [];% need to specify empty inequality constraints

cons = [-hold_pc]; % minus so that pu>=baru

end

%utility constraints fixed product
function [cons, conseq] = constr2FP(x,baru, beta_i,alpha_i,mProdF)

%reshape product
mX_c=mProdF'; % one charac per column
vr_c=x;

Mutil=util(beta_i,alpha_i,mX_c,vr_c) ;% get utilities
pu = PU(Mutil);% get expected utility

hold_pc= pu-baru; %need to be a row vector

conseq = [];% need to specify empty inequality constraints

cons = [-hold_pc]; % minus so that pu>=baru

end


%generate the survival probabilities
function theta_ic = theta(vIncome, vr_c, beta_I, beta_r)
sI=size(vIncome); % vIncome need to ne a row vector
sr=size(vr_c); % vr_c need to ne a row vector
d_i=repmat(beta_I*vIncome,1,sr(1));
d_c=repmat((beta_r*vr_c)',sI(1),1);
theta_ic=1-(d_c+d_i);
end


%generate the loan size
function L_ic = Loan(vIncome,phi, alphatilde, vbetatilde, mX_c,vr_c)
sI=size(vIncome); % vIncome need to ne a row vector
sr=size(vr_c); % vr_c need to ne a row vector
L_i=repmat(phi*vIncome,1,sr(1));
L_c=repmat(((vbetatilde'*mX_c)'+alphatilde*vr_c)',sr(1),1); %vbetatilde is a row vector mX_c contain a product per columns
L_ic=exp(L_c'+L_i);
end


%generate the loan size for FP
function L_ic = LoanFP(vIncome,phi, alphatilde, vbetatilde, mX_c,vr_c)
sI=size(vIncome); % vIncome need to ne a row vector
sr=size(vr_c); % vr_c need to ne a row vector
L_i=repmat(phi*vIncome,1,sr(1));
L_c=repmat(((vbetatilde'*mX_c)'+alphatilde*vr_c)',sI(1),1); %vbetatilde is a row vector mX_c contain a product per columns
L_ic=exp(L_c+L_i);
end



%generate the marginal costs
function mc_c = MC(vmc,vmc2, bankmc, mX_c)
mc_c=(bankmc+ vmc'*mX_c+vmc2' *mX_c.^2)';
end

%generate the utilities
function u_ic = util(beta_i,alpha_i,mX_c,vr_c)
u_ic=beta_i*mX_c+alpha_i*vr_c';%matrix with preference of i in each row product c in each column
end

%generate the market shares
function ms_ic = marketshare(Mutil)
su=size(Mutil);   % to get the number of borrowers
expU=exp(Mutil);
normalization=max(expU')'; %get the maximum of each row
expU=expU./repmat(normalization,1,su(2));
sumexpU=expU*ones(su(2),1);
MsumexpU=repmat(sumexpU,1,su(2));
ms_ic=expU./MsumexpU; %matrix vith preference of i in each row
end


%generate the promised utility constraint
function pu = PU(Mutil)
su=size(Mutil);             % to get the number of borrowers
expU=exp(Mutil);
normalization=max(expU')'; %get the maximum of each row
expU=expU./repmat(normalization,1,su(2));
pu=expU*ones(su(2),1);
end


